#set working directory, clear workspace
  rm(list=ls())

#installing/loading necessary pacakges
#install.packages("R2WinBUGS")
#install.packages("foreign")
#install.packages("mcmcplots")
#install.packages("plyr")
  library(coda)
  library(foreign)
  library(R2WinBUGS)
  library(mcmcplots)
  library(plyr)

##################################################################################################################################  

#reading in matrix of votes and economic variables
  dat <- read.csv("hungary-vote-data-aggregate.csv")

#defining the inflation gap as difference between two year forecast and current target, deleting observations with NA
  infgap <- dat$twoyearfc-dat$target
  dat <- cbind(dat,infgap)
  dat <- na.omit(dat)

#defining variable names to match BUGS file
  y <- dat$currentrate+dat$outcome
  nvote <- length(y)
  #bank <- dat$cbid
  #nbank <- max(bank)
  sig <- dat$vote.unc
  lag <- dat$currentrate
  ygap <- dat$ygap
  igap <- dat$infgap
  xr <- dat$xr

#reading in the model file
  model.file <- "rcm-corr-agg.txt"

#setting intial values
  inits <- function() { list(a=rnorm(1),
                           bsig=rnorm(1),
                           blag=rnorm(1),
                           bygap=rnorm(1),
                           bigap=rnorm(1),
                           bxr=rnorm(1),
                           
                           sigma.a=runif(1),
                           sigma.bsig=runif(1),
                           sigma.blag=runif(1),
                           sigma.bygap=runif(1),
                           sigma.bigap=runif(1),
                           sigma.bxr=runif(1))  }

#running model in WinBUGS
  bugs.model <- bugs(data=c("y","nvote","sig","lag","ygap","igap","xr"),
                    inits=inits,
                    parameters.to.save=c("a","bsig","blag","bygap","bigap","bxr",
                                         "mu.a","mu.bsig","mu.blag","mu.bygap","mu.bigap","mu.bxr",
                                         "sigma.a","sigma.bsig","sigma.blag","sigma.bygap","sigma.bigap","sigma.bxr"),
                    model.file= model.file,
                    n.chains=3,
                    n.iter=900000,
                    n.burnin=300000,
                    n.thin=5,
                    debug=TRUE)

#Extracting Coefficients from bugs.model for Table 1:
  bugs.model #Deviance, DIC
  bugs.model$mean$mu.a
  bugs.model$sd$mu.a
  bugs.model$mean$mu.bsig
  bugs.model$sd$mu.bsig
  bugs.model$mean$mu.blag
  bugs.model$sd$mu.blag
  bugs.model$mean$bygap
  bugs.model$sd$bygap
  bugs.model$mean$bigap
  bugs.model$sd$bigap
  bugs.model$mean$bxr
  bugs.model$sd$bxr
  
  caterplot(mcmcout=bugs.model,"bsig",quantiles=list(outer=c(0.025,.975),inner=c(0.05,.95)))
  bugs.mcmc <- as.mcmc(bugs.model)
  cred.int95 <- HPDinterval(bugs.mcmc,prob=.95)
  cred.int99 <- HPDinterval(bugs.mcmc,prob=.99)
